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Abstract 

The properties of systems with Bose-Einstein condensate in external time-independent 
random potentials are investigated in the frame of a self-consistent stochastic mean- 
field approximation. General considerations are presented, which are valid for finite 
temperatures, arbitrary strengths of the interaction potential, and for arbitrarily strong 
disorder potentials. The special case of a spatially uncorrelated random field is then 
treated in more detail. It is shown that the system consists of three components, 
condensed particles, uncondensed particles and a glassy density fraction, but that the 
pure Bose glass phase with only a glassy density does not appear. The theory predicts 
a first-order phase transition for increasing disorder parameter, where the condensate 
fraction and the superfluid fraction simultaneously jump to zero. The influence of 
disorder on the ground-state energy, the stability conditions, the compressibility, the 
structure factor, and the sound velocity are analyzed. The uniform ideal condensed 
gas is shown to be always stochastically unstable, in the sense that an infinitesimally 
weak disorder destroys the Bose-Einstein condensate, returning the system to the nor- 
mal state. But the uniform Bose-condensed system with finite repulsive interactions 
becomes stochastically stable and exists in a finite interval of the disorder parameter. 
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1 Introduction 



The existence of the condensate fraction and its relation to the superfluid fraction in random 
Bose media have been an intriguing research subject for many years. First, the objects of 
interest have been ^He-filled porous media, such as Vycor glasses, aerogel glasses, and grained 
powders [1,2]. Recently, the physics of dilute Bose gases has gained much interest (see the 
books [3] and review articles [4-7]). Several experiments with Bose- Einstein condensates 
in random potentials have been accomplished, and different techniques of creating random 
fields have been proposed. For example, random potentials can be formed by laser speckles 
[8,9] or by randomly-varying magnetic fields in the close proximity of a current-carrying 
wire [10]. Quasi-random potentials can also be created by using two-color quasiperiodic 
noncommensurate optical lattices [11]. 

In the theory of disordered Bose systems, one considers two types of models. Of one 
type are the lattice models characterized by a boson Hubbard Hamiltonian with random site 
potentials. Such random potentials suppress or may even can completely destroy the long- 
range order related to Bose- Einstein condensates [12]. Fisher et al. [13] have suggested that 
sufficiently strong disorder in a lattice leads to the appearance of a new phase, different from 
insulating and superfluid phases. This is the Bose glass phase, which is characterized by a 
finite compressibility, the absence of a gap in the single particle spectrum, and a nonvanishing 
density of states at zero energy. The phases in these lattice models can be classified [14] on 
the basis of two order parameters, the condensate fraction no and the superfluid fraction Us- 
In the insulating phase, tiq = and Us = 0. For the superfluid phase, both order parameters 
are nonzero, uq ^ and Us 7^ 0. And for the Bose glass phase, there is uq 7^ 0, but there 
is no superfluidity, rig = 0. The occurrence of the lattice Bose glass, arising between the 
insulating and superfluid phases, has been investigated in several theoretical papers [15-21] 
and confirmed in a recent experiment [11]. 

In a second class of models the disordered bosons can be thought of as being immersed in 
an initially uniform system in a random external potential, with no regular lattices imposed. 
This type of models was first studied by Huang and Meng [22], who considered the case 
of asymptotically weak interactions and of asymptotically weak disorder in the Bogolubov 
approximation. Their results were recovered by Giorgini et al. [23] using the hydrody- 
namic approximation, which is mathematically equivalent to the Bogolubov approximation. 
Lopatin and Vinokur [24] estimated the shift of the critical temperature due to weak disor- 
der in a weakly interacting gas, which also was studied by Zobay [25] , using renormalization 
group techniques. If the results obtained for asymptotically weak disorder are formally ex- 
tended to strong disorder, then one comes [22,26] to the state, where no 7^ but = 0, 
which corresponds to the Bose glass phase. However, Monte Carlo simulations [27] for a 
gas with strong disorder, although it confirmed that the superfiuid fraction can be smaller 
than the condensate fraction, found no presence of the Bose glass phase. Also, no Bose 
glass was found in the random-phase approximation at zero temperature and asymptotically 
weak interactions [28]. Instead, increasing disorder led to a first-order transition from the 
superfluid to the normal phase. Thus, the situation with Bose-condensed systems in random 
potentials is well understood for the limit of weak interactions and weak disorder. However, 
it remains controversial when the interactions and/or the disorder become larger. 

The aim of the present paper is to develop a new approach for treating Bose-condensed 
systems in random potentials, when particle interactions and strength of disorder can be 
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arbitrary. We analyze the main properties of the system and the influence of disorder and 
the interaction strength on these properties. In particular, the ideal uniform gas with Bose- 
Einstein condensate is shown to be stochastically unstable, in the sense that an inflnitesimally 
weak random noise destroys the condensate, turning the system to the normal noncondensed 
state. The stochastic instability could be one of the reasons why the ideal Bose-Einstein con- 
densation is not experimentally possible, and confining potentials and atomic interactions 
are necessary for the Bose-Einstein condensation to be realized in the laboratory. Nonvan- 
ishing repulsive atomic interactions stabilize the condensate, which can then exist in a finite 
domain of temperatures and of the disorder strength. At a temperature-dependent value of 
the latter the Bose-condensed system undergoes a first-order phase transition and transforms 
to the normal phase. 

Throughout the paper a system of units is used, where h — 1 and — 1- 

2 System Hamiltonian 

The Hamiltonian energy operator is taken in the standard form 

^[^] = I ^t(r) - ^ + e(r) V^(r) + ^ $o / V^t(r)^t(r)^(r)^(r) dr , (1) 

in which '0(r) = ■^'(r, t) is the Bose field operator, ^(r) is a random external potential, and 
the particle interaction strength 

$0 = 47r ^ (2) 

m 

is expressed through the scattering length and particle mass m. 

The averaging over the random potentials will be denoted by the double angle brackets 
^ . . . The distribution over the random fields is assumed to be zero-centered, so that 

« e(r) » = . (3) 

The stochastic average 

« a^ W) » = ^(r - r') (4) 

defines the correlation function R{r) . The random potential and the correlation function are 
supposed to be real and the latter is also symmetric, such that 

r(r) = e(r), R%r) ^ R{-r) ^ R{r) . (5) 

Therefore their Fourier transforms enjoy the properties 

Cfc = C-fc , R*k = R-k = Rk ■ (6) 

The Fourier transform possesses also the important property ■Cfc ^ Oi when — > oo as 
explained in the Appendix A. In the Fourier representation, Eq. (4) reduces to 

« Ck^p » = hpRkV . (7) 

For the particular case of white noise, when 

R{r) = Ro6{r) , (8) 
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one has 

« CkCp » = SkpRoV . (9) 

The main part of the present paper will not depend on the particular type of the dis- 
tribution over the random potentials, and hence on the concrete choice of the correlation 
functions (4) and (7). But at the final stage, in order to illustrate practical calculations, we 
shall specialize to the white noise characterized by Eqs. (8) and (9). 

All operators from the algebra of local observables are functionals of the field operators 
^(r) and t/'^(r) and of the random variable C,{r). This implies that there are two kinds of 
averages. One kind is the stochastic average <^ A ^ over the distribution of the random 
potentials. And another one is the quantum average with respect to a Hamiltonian H, which 
is denoted as 

<A>H= Trpi, (10) 

with the statistical operator 

Here the Hamiltonian H includes, but is, in general, different from H and remains to be 
specified below. /? = 1/T is the inverse temperature and the trace is over the Fock space 
generated by the related field operators [29,30]. The total average will be denoted as 

< i > = < Tr pi > . (12) 

To describe a Bose-condensed system, where the global gauge symmetry is broken, one 
employs the Bogolubov shift 

^(r) ^(r)=r/(r)+^i(r), (13) 

where ?7(r) is the condensate wave function. The field variable ^^(r) and the operator ipi{r) 
are taken as linearly independent and orthogonal to each other, 

J 77*(r)V'i(r) dr^O. (14) 

ipi{r) is the operator of uncondensed particles, satisfying the Bose commutation relations 
[31-33]. The condensate function is normalized to a fixed, still undetermined, positive value 
Nq, the number of condensed particles 



J\v{r)\'dr. (15) 



The physical value of Nq must then be chosen by minimizing the thermodynamic potential. 
The number of uncondensed particles Ni = N — Nq is given by the average 

Ni = < TVi > (16) 

of the number-of-particle operator 

N, = I Vl(r)Vi(r) dr . (17) 



4 



The total number of particles in the system is 

N = <Ar> = No + Ni, (18) 

with the operator 

N = J Vi^(r)Vi(r) dr ^ No + Ni , (19) 

in which ■ip{r) is the shifted field operator (13). 

According to these definitions, for the correct description of a Bose-condensed system, 
which would be self-consistent in any approximation, one therefore has to employ a repre- 
sentative ensemble [34] taking into account the normalization conditions (15) and (16) or 
(18). This requires [34-36] to use the grand Hamiltonian 

H = H- noNo - iiiNi , (20) 

where H = Hl^j], while fiQ and are the Lagrange multipliers guaranteeing the validity 
of normalizations (15) and (16). Here we shall consider an equilibrium system, but a simi- 
lar representative ensemble can also be defined for nonequilibrium Bose-condensed systems 
[34,37]. 

3 Thermodynamic Potential 

For the frozen disorder, the grand thermodynamic potential is 

Q = -T < InTre-^^ > . (21) 

To provide thermodynamic stability, potential (21) is to be minimal with respect to the 
number of condensed particles. 

The system free energy can be defined as 

F = Q + iioNo + fiiNi . (23) 

At the same time, keeping in mind that in standard experiments only the total number of 
particles N is fixed, but not Nq and separately, we may write 

F^n + I^N . (24) 

Comparing Eqs. (23) and (24) yields the definition of the system chemical potential 

H = hqUq + niui , (25) 

in which no = Nq/N, rii = Ni/N are the corresponding fractions of particles, satisfying the 
normalization condition no + ni — 1. 

It is worth noting that, instead of working with the grand ensemble containing two 
Lagrange multipliers, we could resort to the canonical ensemble with no Lagrange multipliers 
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but with two constraints that are to be satisfied at each step of any calculational procedure. 
One constraint is that the number of condensed particles A^o — ^o{T, N) be fixed by stabihty 
conditions, while the total number of particles N be kept fixed at each step, but not solely 
on average. Such a canonical ensemble could, probably, be realized with the help of the 
Girardeau- Arnowitt representation [38]. However, a weak point of the latter is not only 
that it leads to rather cumbersome calculations but, most importantly, that it does not 
allow simple self-consistent approximations. For instance, it is well known that the Hartree- 
Fock-Bogolubov (HFB) approximation is not self-consistent in the frame of the Girardeau- 
Arnowitt representation, yielding an unphysical gap in the spectrum [38] for a uniform Bose 
system. Girardeau [39] stressed the necessity to deal with the complete Hamiltonian in order 
to make the canonical-ensemble approach self-consistent and to remove the unphysical gap. 
Indeed, Takano showed [40] that this could really be done at least in principle, if one would 
use all terms of the Hamiltonian. However this necessity makes the problem practically 
unsolvable: in general, an exact solution for the problem is not known, and as soon as an 
approximation is involved, one confronts the danger of getting not self-consistent results 
[41]. Contrary to this, relaxing the imposed constraints, by introducing the corresponding 
Lagrange multipliers, being mathematically equivalent, makes all calculations much simpler, 
at the same time preserving the theory self-consistency for any given approximation [34-37] . 

In order to calculate the thermodynamic potential (21) for the frozen disorder, one often 
takes recourcc to the so-called replica trick, as is used in the theory of spin glasses [42]. Here 
we shall employ another approach, based on the method of separation of variables. The idea 
of this method is as follows. The main aim is to transform the given Hamiltonian if to a 
separable form 

Hsep — Hq + , (26) 

in which Hg depends only on quantum variables, while depends only on classical stochastic 
variables. Such a transformation can be achieved by means of canonical transformations and 
some simplifications. Then the corresponding thermodynamic potential 

Qsep = -T < In Tre-^^=«^ > (27) 

reduces to the sum 

aep = -TlnTre-^^«+<i/^> , (28) 

in which the manipulations with quantum and stochastic variables are separated. If the 
separable Hamiltonian (26) does not exactly represent the initial H, so that 

H = H,,p + h {h = H- Hsep) , (29) 

then corrections to the thermodynamic potential can be obtained by perturbation theory 
with respect to h, giving in the second order 

n = nsep + <h> -pA'^Ch) , (30) 

where A'^{h) is the dispersion A^(/i) =< fi^ > — <h . In agreement with Eq. (30), one 
has VL < VLsep+ < h >, which is the Gibbs-Bogolubov inequality. 

The method of separation of variables has no need for the replica trick. The derivation of 
the separable Hamiltonian (26) can be accomplished by means of decouplings and canonical 
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transformations and does not require the existence of small parameters. All essential non- 
linearities with respect to particle interactions and disorder strength can be preserved in the 
Hamiltonian (26). The use of the Gibbs-Bogolubov inequality, mentioned above, can be done 
in the standard variational way, by minimizing the right-hand side of this inequality, which 
again does not require the existence of small parameters. Therefore this method makes it 
possible to consider strong interactions and strong disorder. 



4 Stochastic Quantization 

According to Eq. (3), the external random potential is zero on average. This allows us to 
treat the condensate wave function, which is the system order parameter, as uniform, so 
that [22] 

^(r) =<^^(r)>= (31) 

where i^{r) is the shifted field operator, po = Nq/V is the condensate density, and the total 
average (12) is assumed. In agreement with Eq. (13), one has 

< V'i(r) > = . (32) 

Expanding the field operators of uncondensed particles in plane waves, we represent the 
grand Hamiltonian (20) as the sum 

/^=E^^"^+^exf (33) 
n=0 

Here the zero-order term 

//W = Q po^o - /^o) No (34) 

does not contain the operators of uncondensed particles. For the first-order term, because 
of the property (14), we get 

= . (35) 
The term of second order, with respect to the operators a^, becomes 

^^'^ = E + 2po^o - /^i j alak + ^ E Po^o (ajal^ + a-^ak) ■ (36) 
For the third-order term, we have 

^(3) ^ ^ (alak+pa-p + alpal^pak) , (37) 

k,p 

where the prime on the summation symbol implies that k 7^ 0, p 7^ 0, k -|- p 7^ 0. The 
fourth-order term is 

H^'^ = 4? E E' $o4aJa,_,a,+, , (38) 

Q k,p 
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where the prime on the summation sign means that k7^0,p7^0,k — q7^0,p + q7^0. 
The last term in Eq. (33) corresponds to the action of the external random field, given by 
the expression 

Hext = Poio + i^l^k + Ck(^k) + ^ Yl o.Wik-p ■ (39) 

When one assumes asymptotically weak interactions, one omits the terms /J'^^^ and H'^'^\ 
thus, coming to the Bogolubov approximation [31-33] . Since we aim at considering arbitrar- 
ily strong interactions, we have to keep all terms of Hamiltonian (33). But we may simplify 
the terms H^^^ and H^^^ by means of the Hartree-Fock-Bogolubov (HFB) approximation 
[35,36]. Then, we get 

H^^^ = . (40) 
To express the result for the term H^^^ in a compact form, we introduce the normal average 

Uk =< a\ak > , (41) 

which is the momentum distribution of atoms, and the anomalous average 

=< o-ka-k > ■ (42) 

The quantity |o"fc| can be interpreted as the momentum distribution of paired particles [35]. 
Then the density of uncondensed particles is 

Pi = Y.^k, (43) 

while the sum 

= TT E ^fe (44) 

gives the density \ai\ of paired particles. Applying the mean-field approximation we find 
from Eq. (38) 



+ — E "^0 nk+pa^ap + - (ak+palaip + al^pU^pUp^ - - (rik+pHp + (Jk+pCr*^ 



(45) 



A special care has to be taken in reorganizing expression (39) describing the interaction 
of atoms with external random fields. The second term in Eq. (39) corresponds to linear 
interactions between random fields and atoms, while the third term describes nonlinear 
interactions. If one omits the third term, as has been done by Huang and Meng [22], thus, 
keeping solely the linear interactions, then one limits oneself by weak disorder. Since our 
aim is to consider arbitrarily strong disorder, wc need to keep this term. The difficulty with 
treating the nonlinear term in Eq. (39) is that, in the mean- field approximation, it is zero 
on the average, as far as 

< alap^k-p > ~ < alap >< ^k-p > = . (46) 
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If we would treat this term in the simple mean-field manner replacing alap^k-p by 

< al.ap > ^k-p + OfeOp < Cfe-p > — < Oj^Op >< ^k-p > — ^kp^k^o , 

we would kill all quantum effects, reducing the term to the trivial form. The way out of this 
problem is to employ a more refined approximation. 

We shall use the ideas of the stochastic m,ean-field approximation, which has been applied 
to accurately treat quantum and stochastic effects in systems interacting with electromag- 
netic fields [43] and in spin systems [44,45]. In considering these systems, one encounters 
the same type of the difficulty. If one uses the simple mean-field approximation, often 
called semiclassical, then quantum and random effects are washed out, which may lead to 
principally wrong results. To accurately take account of the latter effects, the mean-field 
approximation is to be modified [43-45]. 

Let us remember that we have two types of averages for any operator A. The stochastic 
average <^ and the quantum average < A >h defined in Eq. (10). The operators of 
uncondensed particles and aj. are, strictly speaking, functions of the random fields ^k- We 
may separate the quantum and stochastic averages and consider the quantum average 

ak =< ak>H , (47) 

which is a function of the random fields. This quantity ak is not zero, even though its total 
average 

<afc> = <afe> = (48) 

is of course zero, according to Eq. (32). In the nonlinear term of Eq. (39), in the spirit of 
the stochastic mean-field approximation [43-45], we now make a mean-field type decoupling 
with respect to the quantum averaging only, not with respect to the stochastic average, that 
is, we write 

alttp^k-p = (4 < ap>H + < al >h ap- < aj. >h< ap >h) Cfe-p ■ (49) 

One may notice that if we would employ in decoupling (49) the total averages of type (12), 

instead of the quantum averages of type (10), then the left-hand side of Eq. (49), according 
to Eq. (48), would be reduced to zero, similar to Eq. (46). In order to retain the influence 
of the left-hand side term of Eq. (49), we invoke here not the total but only the quantum 
averages. Using the latter, instead of the total averages (12), makes decoupling (49) more 
general, thus, allowing us to retain the influence of nonlinear stochastic terms [43-45]. Let 
us also deflne the stochastic held 



<fk = ^k + jrYl ^pik-p ■ (50) 

^ ^ P^O 



Then the random-fleld Hamiltonian (39) transforms to 

Hext = PoCo + Yl i^k'Pk + 'Plak) - ^ H ttfettpCfc-p , (51) 



k^o ' ^ k,p{¥^o) 



where Co = / C(r)'ir. 
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Finally, introducing the notation 



Uk = — + 2p$o - A^i , (52) 

where p = po + Pi is the total particle density, and defining 

A= (po + (7i)$o , (53) 
we obtain for Hamiltonian (33) the form 

H = Ehfb + X] ^fe4«ifc + ^ H ^ + O-feOfe) + Hext , (54) 
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in which 



- ^ (2p? + a?) ^oN (55) 



Ehfb 

2p 



and i^ext is given by Eq. (51). It is worth emphasizing that Hamiltonian (54) has the mean- 
field form with respect to the field operators a^, but it contains, via Hg^^, the nonlinear 
terms with respect to the random variables ^k, CKk, a-nd (fk- The latter allows us to consider 
disorder of arbitrary strength. 



5 Separation of Variables 

Quantum and stochastic variables in the Hamiltonian (54) are yet intermixed. To separate 
them, we shall use the method of canonical transformations. First, we employ the usual 
Bogolubov canonical transformation 

Ofe = Ukh + v*_kb-k ' «fe = K^l + v-kb-k ■ (56) 
Using these in Eq. (54), we get 

H = EB + ^Skbibk + H,,t: (57) 

where 

Eb = Ehfb + ::J2i^k - ouk) (58) 
and £jt is the Bogolubov spectrum 



Sk = ^/^4^ ■ (59) 
Equation (51), containing random fields, now becomes 

Hext = PoCo + Yl i^k^k + Dlhk) - ^k<^pik-p , (60) 

where 

Dk = {ul + vl)ipk. (61) 
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The coefficient functions in transformation (56) are defined by the equations 

Then we apply another canonical transformation 

bk^h - bl^bl- ^, (63) 

which transforms Hamiltonian (60) into 

H^EB + Y.^kblh + H^, (64) 

where Eb is the nonoperator part (58), the second term does not depend on stochastic 
variables, while the last term 

contains only stochastic fields, but no quantum variables. 

In that way, the quantum operator variables 6^ and b], and the stochastic fields ^k, Oik, 
and ifk are separated in Hamiltonian (64). This will allow us to calculate different averages 
and to analyze the influence of random fields on the system. 



6 Random Fields 

Let us, first, consider the Bogolubov spectrum (59). As is seen, it does not explicitly depend 
on the random fields, thus, representing the spectrum of collective excitations for a system 
that is uniform on the average. For a uniform system, there exists the Hugenholtz-Pines 
theorem [46,47] requiring that the spectrum be gapless, so that 

lim£jt = 0, £jt > . (66) 
fc— >o 

Then, from Eqs. (52), (53), and (59), it follows that 

Ati = (p + Pi -(Ji)^o ■ (67) 

As a result, Eq. (52) reduces to 

uj, = ^ + A. (68) 
zm 

The Bogolubov spectrum (59) acquires the form 



.(ckr^{£^\ (69) 
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in which the sound velocity 

(70) 

is expressed through the quantity 

A = mc^ = (po + (Ti)$o , (71) 

following from Eq. (53). 

Another way of deriving Eq. (67) and, respectively, the Bogolubov spectrum (69) 
is as follows. We may consider the equations of motion for the matrix Green function 
G{k, uj) = [Gaj3(k, a;)] as has been done by Bogolubov [48]. The presence of the random-field 
Hamiltonian (39) contributes to these equations with the terms all of which, in the mean- 
field approximation, can be set zero, in agreement with Eq. (46). For the Green functions, 
one has the Bogolubov theorem [48] 

|Gn(k,0)-Gi2(k,0)|>^, 

from which the Hugenholtz- Pines relation fii — Eii(0, 0) — Ei2(0, 0) follows, where S(k, lo) — 
[Ea/3(k, a;)] is the matrix self-energy. In the HFB approximation, we have Eii(0, 0) = 2p$o 
and Ei2(0, 0) = (po + <^i)^o- This gives us exactly the same equation (67). 
Gombining the canonical transformations (56) and (63), we get 

dk = Ukbk + t'-fcO-fc , %= u^bl + v-kb-k ^k ■ ( '2) 

Because of the form of the Hamiltonian (64), one has < bk >—< bkbp >— 0. Then, from 
Eqs. (47) and (72), we find 

^k f„n\ 
Wfc + A 



Hence, 



« » = , ' , . (74) 



|2 

{uk + A)2 

By Eq. (48), we also have < <pfc >—<^ ^pk 0. Substituting relation (73) into Eq. (50), 
we come to the equation 

^. = ^a-^E%^ (75) 
V a;p + A 

defining the random field ip^- This is a Fredholm equation of the second kind. 

Using Hamiltonian (64), it is straightforward to get the momentum distribution of quasi- 
particles 

TT, ^ < blk > = (e^^'= - l)"' = icoth (^) - ^ ■ (76) 
For the momentum distribution of atoms (41), we find 

Uk = (ul + vl) TTk + vl+ <^ \akf :^ (77) 
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and for the anomalous average (42), we have 



CTfc = (1 + 2nk)ukVk+ < \akf > . (78) 
With Eqs. (62) and (76), we finally obtain the normal average 

= ^ (^) I"""!' ^ ^^^^ 

and the anomalous average 

The contribution of the random potential comes through the last terms in Eqs. (79) and 
(80). These terms are related to the random field ipk by means of Eqs. (73) and (74). And 
the random field (pk is defined as the solution of the Fredholm equation (75). 



7 Glassy Fraction 

In order to elucidate the physical meaning of the terms, induced by the random potential, 
let us draw some analogies with the theory of spin glasses [42]. For the Bose system, we may 
define an order parameter, which is the analogue of the Edwards- Anderson order parameter 
in spin glasses [42]. To this end, we recall that the total average < ■^i >= 0, according to Eq. 
(32). But, separating the quantum and stochastic averages, we can introduce the density of 
the glassy fraction 

= :^y<|<V'i(r)>^|2»(ir. (81) 



Pg 



Passing to the Fourier transform of ■0i(r) and using Eqs. (72), we reduce Eq. (81) to 

1 

V 



PG-hl2^ » . (82) 



Consequently, the meaning of the quantity 



is the momentum distribution of the glassy fraction. We may assume that the nominator 
of Eq. (83) is not increasing with k. However, its denominator, according to Eq. (68), 
increases with k as k'^. Hence, distribution (83) is a rapidly decreasing function of k, with 
its maximum at k — 0, where 

= — ^ . 84 

The glassy density (82), using relation (74), can be represented as 
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Since the integrand in Eq. (85) falls off rapidly and <S \^k\^ ^ is slowly varying with /c, we 
may substitute <^ |</7o|^ ^ instead of <S |</7fcp which gives 

PG = ^ ^^(0) . (86) 

TT 

For the dimensionless glassy fraction, we then have 

PG {mcf 

nG = — = noiO) . (87) 

P Trp 

Let us consider the glassy density matrix 

r dk 
pair,, r2) = j nG(k)e^'^--- ^ , (88) 

in which ri2 = ri — r2. This, with the glassy distribution (83), gives 

Taking into account that the main contribution to integral (89) comes from small k, and 
using the equality 

1"^ a;sin(ax) na 
io (62 + a;2)2 = 46 ' 

we obtain the glassy density matrix 

pair, 0) = pGe-''°' {ko = 2mc) . (90) 

This demonstrates that the locahzed short-range order of the glassy fraction has the decay 
length l/ko, which coincides with the healing length. 

It is important to stress that the presence of the glassy fraction in the type of systems 
under consideration here does not turn the whole system into a Bose glass. This is because 
by the commonly accepted classification, the Bose glass phase requires that the superfiuid 
fraction Ug be zero, which is not the case here. Also, the density of states 

47rk\u) dk{uj) 



p{uj) = 



(27r)3 duj 



in which k{ijj) is defined by the equation = lv, with Sk from Eq. (59) or cquivalently Eq. 
(69), yields 

pH = ^ — ^ . 

This tends to zero at small u; as 

Thus, the system does not represent a Bose glass, for which p(0) must be finite. 

To conclude, the action of external random fields on the Bose system induces the appear- 
ance in the latter of the glassy fraction but need not transform the system as a whole into 
the Bose glass phase. 
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8 Thermodynamic Stability 



It is interesting to study the influence of random potentials on the thermodynamic stability 
condition (22). For the Lagrange parameter /^o of the condensate fraction, introduced in Eq. 
(20), we have from the first of Eqs. (22) 

//o = (p + Pi + (7i)$o + A*G , (91) 

The last term ^ 

'"^ = 9 /ATT/ E < 4a + Ck^^k > (92) 
2ViVo V 

is caused by the direct action of the random potential. Prom the second of Eqs. (22), we 
find 

*o > ^ . (93) 
2po 

Thus, the stabihty condition (93) for the particle interaction strength $o of Eq. (2) depends 

on the value fic- 

Equations (72) and (73) show that 

< al^k > = < a^^k > = ,A ■ 

Thus, the glassy term (92) takes the form 

^ _ 1 < + CkVk > .g^V 

2^ V ^0 oOk + A ■ ^ ' 

One has to exercise considerable caution when analyzing Eq. (94). To stress this, let us 
start with the attempt of calculating fie by means of perturbation theory with respect to 
weak disorder. Under asymptotically weak disorder, the limiting approximate solution of 
Eq. (75) is 

^k - ■ 

Substituting this into Eq. (94) yields the perturbative expression 
With the definition of Rk in (7), we get 

Rk dk. 



-/ 



uJk + A (27r)3 ' 



which is exactly the form obtained in Ref. [24]. Since the correlation function R^ is assumed 
to be positive, one has < 0. Then condition (93) tells us that the action of the random 
potential stabilizes the system, which does not appear plausible, physically, however. 

On the other hand, if one interprets the random potential as being caused by the presence 
of randomly distributed impurities, which then justifies the use of analytic regularization 
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procedures for physical integrals, and if one takes the limit of the uncorrelated spatial white 
noise potential, defined in Eq. (9), then one gets from Eq. (94) the different perturbative 
value 

2 

fj,Q — — m^cRo ■ 

TT 

The latter is positive, contrary to /Iq < 0. In this way, the sign of the glassy term (94), in a 
perturbative evaluation, is not independent of the method of calculation. In other words, it 
remains unclear whether the random potential stabilizes or rather destabilizes the system. 
This gives a strong hint that the application of perturbation theory with respect to weak 
disorder may not be justified for the considered case. This would show up via inconsistencies, 
such as divergencies, when going to higher order in the perturbative calculations we sketched 
here. 

Fortunately, we are able to calculate Eq. (94) without resorting to the weak-disorder 
approximation, but by considering instead the whole Eq. (75) exactly. We immediately 
obtain then 

«y^^ + y^o» «a + eo» . .... 

- 2v/]Vo W ^ ° ■ ^^^^ 

Thus, we find fic = for any type of the random potential and any strength of disorder. So, 
the stability condition (93) acquires the simple form $o > 0. 

This result teaches us that the action of random potentials on Bose systems may lead 
to nonperturbative effects, when calculations for asymptotically weak disorder can yield 
incorrect conclusions. 



9 Energy Contribution 

The direct contribution of the random fields to the internal energy of the system is given by 
the average of term (65) entering the Hamiltonian (64), that is, by 

= < //^ > = < i/^ > . (96) 

With relation (73), the latter gives 

^ ,^0 ^k + A V ^^^^^ {uj, + A){uj, + A)- 

Exercising now the required caution when dealing with random fields, we shall not use 

perturbation theory for weak disorder, but shall instead take into account the exact Eq. 
(75) and use the properties — — in line with Eqs. (6). Employing once Eq. 

(75), we have 

77 < <^k(^pik-p > = —r. — r — ■ 

This allows us to transform Eq. (97) into 

(98) 
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Invoking once more Eq. (75) in the form 



E^pr-p _ V ^ 'U ^ 



we reduce Eq. (98) to 
which results in 



= JNo < (/?o > -Po < Co > , (99) 



E^^O (100) 



for any kind of the random potentials and any strength of disorder. 

It is instructive to stress again that the usage of perturbation theory with respect to 
weak disorder is not appropriate here. Really, if we substitute the approximate solution 



iff, ~ y Po/V^k: corresponding to weak disorder, into Eq. (98), we get the perturbative 
energy 

^ ya;, + A(27r)3 

in the same form as has been obtained by all other authors using the weak-disorder limit. 
This result would seem to tell us that the presence of random potentials diminishes the 
internal energy. 

However, if we interpret the presence of the random potential as the existence of randomly 
distributed scatterers, use the analytic regularization of integrals, and treat the case of white 
noise, then we find 

2 

E'l = - Nm^cnoRo ■ 

^ 71 

Hence, the internal energy would now seem to increase with Rq. However, both mutually 
conflicting perturbative results are at variance with the exact value (100), which is always 
zero. 

Again, as in the previous Sec. VIII, we come to the conclusion that perturbation theory 
with respect to weak disorder can lead to incorrect results. 



10 Uncondensed Particles 

The properties of uncondensed particles arc characterized, first of all, by their density pi of 
Eq. (43) and the anomalous average cti of Eq. (44). Using Eq. (79), the density pi can be 
presented as the sum 

r dk , . 

Pi^ J rik — = Pn + Pg (101) 



of the normal density 
and of the glassy density 



dk 



(27r) 



(102) 



dk 

PC ^ I ncik) ^ , (103) 
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which can be written as in Eq. (85). The normal density (102) can be represented as 



1/2 



/ mc 
coth^x 



dx 



(104) 



With the help of Eq. (80), the anomalous average (Ti, whose absolute value is the density 
of paired particles, can similarly be written as the sum 



dk 



ctn + Pg 



of two terms. The first term is 

CTjv 



--[ — coth ( 



A ,.(ek\ dk 

coth — - - — -, , 
Ek \2Tj (27r)3 ' 



(105) 



(106) 



while the second term is the same glassy density (103). Equation (106) can be rewritten as 



1 f A 
~ 2J Vk 



coth ( ^ ) - 1 



(2^ 



where 



(7o 



/ 



1 dk 



(107) 



(108) 



2 J Sk (27r)3 • 

The integral in Eq. (108) is ultraviolet divergent. This divergence is well known to be 
unphysical, since it is caused by the usage of the contact interaction potential. A general 

way of treating such integrals is as follows. First, one restricts to asymptotically weak 
coupling and applies the technique of dimensional regularization , which is an accurately 
defined mathematical procedure in that limit [5] . Then one analytically continues the result 
to finite coupling. The dimensional regularization gives 



1 dk 



Bk (27r)= 



2m 



In this way, we find for Eq. (108) 



(mc)^ 



Changing the variables of integration, Eq. (107) can be represented in the form 

(^TT^- 1)1/2 



(mc 



2^2 TT^ Jo 



f 

Jo 



, / mc \ 
coth -— - a; — 1 



2T 



dx . 



(109) 



(110) 



At low temperatures, when T/mc^ -C 1, Eq. (104) gives 



Pn 



{mcf ^ {may ( T 



37r2 



12 



and Eq. (110) yields 



CTjv ^ (To 



{mcY 
12 



T 



mc^ 



;iii) 



(112) 
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In the case of weak interactions, such that m(? jT^ <^ 1, where is the critical temper- 
ature 



T = 

c — 



27r 



m 



P 



C(3/2), 



2/3 



Eqs. (104) and (110) lead to 



and, respectively, to 



T\3/2 (mc 



(JaT ~ (Jo - 



37r2 

m^cT 
27r ■ 



(113) 



(114) 



(115) 



The analysis of the behavior of p^v and (Tjv shows that these quantities are characteristic 
of the Bose system without disorder, while the explicit influence of the random potential is 
contained in the glassy density (103). 



11 Superfluid Fraction 



By a general definition, the superfiuid density is the partial density appearing as a response 
to a velocity boost, 

1 d 

Ps = ^z—rp hm — • < >^ , (116) 

where the average of the system momentum = P + Nmv is calculated with the Hamil- 
tonian 

2' 



H^ = H + J ^t(r) ^_^v. V+^j V'(r) dr 



of the liquid moving with velocity v. 

The dimensionless superfluid fraction can be represented as 

n. = - = l-^, (117) 
where Q is the dissipated heat, having for an equihbrium system the form 

< P2 > 

A detailed derivation of Eqs. (117) and (118) can be found, e.g., in Ref. [4]. 

Passing to the Fourier transforms, we have < P^ >= J2k,p0^ ' P) < ^fc^p where 
■hk = alttk- In the HFB approximation, 

< hkfip > = UkUp + SkpUkil + Uk) + 5-kp(jl . 

Then Eq. (118) assumes the form 



^ 1 /■ fc^ / n ^\ dk , , 
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Taking into account Eqs. (79) and (80), we may represent the dissipated heat (119) as the 
sum 

Q^Qn + Qg (120) 

of two terms. Here the first term 

If dk 

^""'^pl sinh2(£fc/2T) (2^ ^ ^ 

is the heat dissipated by normal uncondensed particles. And the second term 

Qo-^l ^^^^ coth m (122) 
2mpJ £jk(cUfc + A) \2T) (27r)3 ^ ' 

is the heat dissipated by the glassy fraction. 
Equation (121) can be rewritten as 

^ {mcf (VI + - Ifl'^x dx 

4-\/2 n'^mp Jo 

VI + sinh2(TOc%/2r) ' ^ ^ 

At low temperatures, such that T/mc^ <^ 1, we get 

Qn^ ' i^) . (124) 



15mp Vmc^ 

And in the limit of weak interactions, when mc^/Tc <^ 1, we find 



- 



(125) 



^TJ C(3/ 
where (,{■) is a Riemann zeta function. 

12 Sound Velocity 

The sound velocity c enters in the majority of the above expressions. The velocity itself is 
defined through Eq. (71), which can be written as 

mc2 = (p-pi + (7i)$o , (126) 
taking into account that po — p — pi- According to Eqs. (101) and (105), we have 

Pi = Pn + Pg , (Ji = (Jn + Pg ■ (127) 
Therefore Eq. (126) becomes 

mc^ ^ {p - Pn + (Jn)% . (128) 
It is convenient to work with the dimensionless fractions 

un = — , cr = — . (129) 
P P 
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Since uq + rii — 1 and ni = un + nc the normalization 



riQ + riN + riG ^ 1 (130) 

holds true. 

Let us define the gas parameter 

7 = ■ (131) 

and the dimensionless sound velocity 

mc 

Then, taking into consideration the interaction strength (2), equation (128) for the sound 
velocity can be reduced to the dimensionless form 

s2 ^ 47r7(l -niv + (j) . (133) 

At first glance it might seem that the sound velocity, being the solution of Eq. (133), 
docs not depend on the glassy fraction uq induced by the random fields. That fraction is 
defined by Eqs. (85) and (87) which give combined 

However, through normalization (130), no influences the condensate fraction no, and the 
latter enters the anomalous fraction cr, thus, influencing the sound velocity through Eq. 
(133). For example, at zero temperature, according to Eqs. (104) to (112), we have 

= = , ^0 = 1 - ^ - riG (7^ = 0) . (135) 

Increasing disorder increases the glassy fraction nc, so, decreases the condensate fraction 
no, which decreases a. At the same time, the normal fraction un also decreases. Since njv 
and (T enter Eq. (133) with opposite signs, their changes almost compensate each other. 
Numerical calculations show that the sound velocity s as a function of the disorder strength 
slightly decreases with the latter. 

13 Structure Factor 

The structure factor of a random system is deflned as the stochastic average 

S{\i) = < ShO<.) > (136) 

of the frozen factor 

-5i/(k) = ^ / [< n(r)n(r') >h - < n(r) >h< n(r') >h] e'''^^''-''^ drdv' , (137) 
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expressed through the quantum averages, in which n(r) = ^^(r)^(r). Note that Eq. (136), 
in the theory of random systems, is called the connected structure factor. With the Fourier 
transform 

Pfe = y n(r) e-*"- dr , (138) 

Eq. (137) becomes 

ShO^) = ;^ (< PkPk >H - <Pk >H< h >h) ■ (139) 
Invoking the Bogolubov shift (13), for Eq. (138), we have 

Pk = koNo + ^ + sJNq (a^ + Ofc) . (140) 

The quantum averaging of Eq. (140) gives 

<Pk>H^ SkoNo + 1^ < 4+p«P >H + \/No (alfe + Qijfc) , (141) 

where is defined in Eq. (47). Calculating the first term in Eq. (139), we arrange the 
operator product in the normal form and use the second-order procedure, following the 
standard calculations, the same as for Bose systems without disorder [7,35]. Then for the 
structure factor (136), we find 

S{k) = 1 + 2{nk + CTfc) - 4 < \ak\^ > . (142) 

Substituting here Eqs. (79) and (80), we obtain 

m - ^ ^oth . (143) 



2mek \2T 

The central value of the structure factor is known to be related to the isothermal compress- 
ibility 

1 (8V\ 1 (dp\ S(0) 

where P is pressure. To emphasize the role of the glassy fraction, the central structural 
factor can be written as 

S{0) = ^ + Ana, (145) 

where A = — 2TcQ/mcQ; Cq is the sound velocity in a system without disorder, and Cq = 
dc/dnc at the value nc = 0. From numerical calculations it follows that the coefficient of A 
is positive. Thus, the above expressions show that the random field, via inducing the glassy 
fraction no, leads to an increase of the density fluctuations, the isothermal compressibility, 
and the structure factor. The physics of these results seems to be clear. An additional 
external random potential should lead to the increased scattering of either light or neutrons, 
which is characterized by an increase of the structure factor. 



22 



14 White Noise 



The influence of tlie random potential on physical characteristics comes through the corre- 
lator -C l^^fcP S>. To calculate the latter explicitly, we need, first, to solve the random-field 
equation (75) and, second, to specify the type of the random potential, which till now has 
been arbitrary. 

Let us consider Eq. (75) assuming that in the sum of its second term the main contribu- 
tion comes from the region of small momenta (see Appendix A), so that this equation can 
be represented as 

^. = ^&-^E^- ("6) 

This is the Fredholm equation of the second kind with a separable kernel. Such an equation 
can be solved exactly. The corresponding exact solution is 



Calculating <^ \(pk\'^ with Eq. (147), we can use the expansion 

-1 oo 

which requires the knowledge of the stochastic correlators such as -C ^fci^fc2 • • • ^fc„ 

To define these correlators explicitely, we consider the case of the Gaussian white noise 
[49]. Then we obtain 



/ m Ft \ 

«k.r»= Poi?oE(2^ + l)!!(^) ' (148) 

n=0 ^ ^"^^ ^ 

where the integral 



/ 



dk m 



(a;fe + A)2 (27r)3 47rc 

has been used. Note that the right-hand side of Eq. (148) does not depend on k. This allows 
us to find the explicit expression for the glassy density (85), which becomes 



PG = ^ < Wk? > , (149) 

where the right-hand side is given by the series (148). 

As is clear from its form, series (148) is asymptotic with respect to the parameter 
mRo/Airc. In order to define the quantity <^ Iv^fcP ^ for finite values of the latter pa- 
rameter, it is necessary to employ a resummation procedure for series (148). For example, 
we could resort to the Pade summation [50]. Here we shall use another, more general and 
accurate method, based on the self-similar approximation theory [51-53]. We shall make 
use of the method of self-similar factor approximants [54-56]. This method was shown to 
be more general than that of Pade approximants and, contrary to the latter, being uniquely 
defined. The method we use is sketched in the Appendix B. 
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For convenience, we introduce the dimensionless noise parameter 

Then, representing the sum of series (148) by the self-similar factor approximant of second 
order, we obtain for the glassy density (149) the expression 

- ls^l\s - vfn ' ^^^^^ 

in which s is the dimensionless sound velocity (132). 

Taking into account normalization (130), we find the condensate fraction 

754/7 _ ^N3/7 

and the glassy fraction 

^« = . + 7.^A(.-.)3A ' (1^3) 

which are expressed through the normal fraction n^v = Pn/ Pi with p^v given by Eq. (104). 

The case of weak disorder corresponds to a small noise parameter (150). Then the 
condensate fraction (152) is 

"«-('-^-|?)(^-"-' (^^*) 

and the glassy fraction (153) becomes 

^ ( 2z/\ , , , , 

ncc^- 1 + — (1-n^), (155) 



7s V Is) 
when 1/ ^ 1. 

If, in addition, atomic interactions are asymptotically weak, such that — > 0, then the 
glassy fraction (155) tends to 

In this limit, the sound velocity acquires the Bogolubov form 



c ~ = — ^Trpas . (157) 

V m m 



As a result, the glassy fraction (156) transforms to 



pas 



which exactly coincides with the expression found by Huang and Meng [22] in the limit of 
asymptotically weak interactions and weak disorder. 
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We may notice that the noise parameter v enters Eqs. (151) to (155) in the combination 



V 7mRo 



s Attc 



It would, therefore, be tempting to consider the ratio u/s as a new parameter. However, 
this ratio becomes reaUy a parameter solely for asymptotically weak interactions, when 



But at finite interactions, the sound velocity c = c(T, p, a^) is a complicated function of 
temperature, density, and scattering length. Respectively, the dimensionless sound velocity 
s = s{T,p,'~f) is a function of temperature, density, and the gas parameter, defined by Eq. 
(133). Hence, at finite interactions, temperatures, and disorder strength, the situation is 
more involved and one cannot reduce the consideration to dealing with the ratio u/s, which 
is not anymore a parameter. 

Equations (152) and (153) show that when atomic interactions are switched off, so that 
s — > 0, then there are no positive solutions for the fractions no and ric for any finite noise 
parameter i>. This means that the ideal Bose-condensed gas is stochastically unstable, in the 
sense that any infinitesimally weak disorder completely destroys the Bose-Einstein conden- 
sate, rendering the system to the normal state. 

In the case of an interacting Bose-condensed system with a finite gas parameter 7, the 
system is stable below a critical noise parameter Uc = UciT, p,'~f). Increasing z/ diminishes 
the condensate fraction but increases the glassy fraction. Reaching the critical value i^c, the 
system undergoes a first-order phase transition, when uq and uq jump to zero, after which the 
normal phase prevails with rij^ — 1. This is in agreement with a first-order phase transition 
found in the particular case of zero temperature and asymptotically weak interactions [28]. 
According to our numerical estimates, the jumps of no and rig are close to those found in 
Ref. [28] at the transition point i>c- 

When disorder is absent, the system displays the second-order phase transition at the 
critical temperature (113) coinciding with that of the ideal Bose gas, which follows from ex- 
pansions (114) and (125). As soon as there appears disorder, with any finite noise parameter 
z/, the phase transition becomes of first order. At asymptotically small z/ — > 0, numerical 
estimates give the shift of the critical temperature 6Tc ~ — 2i//97r, which is close to the shift 
found in Refs. [24,25]. 

To analyze the behavior of the superfluid fraction (117), we need to know the dissipated 
heat (120). The part of this quantity, due to normal particles, is given by Eq. (123). Another 
part, caused by the heat dispersed by the glassy fraction, is defined by Eq. (122). With the 
white-noise relation (149), expression (122) can be represented as 







{as ^ 0) . 





(159) 



where the notation 




(160) 



is introduced. 
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Integral (160) diverges for any finite £, so that one has to invoke some regularization of 
this integral. There are several ways to regularize the integral, all of which yield the same 
result. 

First of all, we understand that the divergence of the above integral is caused by the 
white noise. For a colored noise, we should go back to Eq. (122), in which <^ \^k\^ ^ would 
be a diminishing function of A;, but not a constant in A;, as for the white noise in Eq. (148). 
Then integral (122) would be convergent. For such a colored noise, we could approximate 
Eq. (122) as 



nfio I 



2mp 



; COth — ; 



dk 

(2^ 



-/ 
mp J 



k"^ < \(Pk\^ > dk 



elicv. + A) (27r)3 



Passing after this to the white-noise relation (149), we obtain 

Qg = 2nGT . 



(161) 



Another way of regularizing Eq. (160) is a kind of analytic regularization, where one, 
first, takes the integral for that region of the parameter e, where the integral converges, and 
then analytically continues the result to arbitrary values of this parameter. The sole domain 
of £, where integral (160) converges, is the region of asymptotically small £ ^ 0. Then 



£ Jo 



°° dx 



TT 

4£ 



(162) 



Substituting this into Eq. (159) gives again Eq. (161). 

We can also apply for integral (160) a resummation regularization, when the integral, first, 
is represented as a series, after which the series is reorganized with the help of a resummation 
procedure. Two variants of the self-similar regularization are described in the Appendix C, 
both of which lead to the same answer (161) as two previous regularizations considered above. 
It is important to stress that wc have accomplished several ways of regularizing integral (160) 
in order to prove that the result does not depend on the regularization procedure involved. 

Thus, combining Eqs. (117), (120), and (161), we find the expression for the superfiuid 
fraction 

* (163) 



1 - 



no - 



3 3T ' 

in which uq is defined in Eq. (153) and Qat, in Eq. (123). It is worth stressing that, though 
the superfiuid fraction (163) is linear with respect to the glassy fraction uq, it is far from 
being linear with respect to the strength of disorder u, as follows from expression (153) for 

the glassy fraction. 

At low temperatures, when T -C mc^, the superfiuid fraction (163) is 



4:U/3 



V + 7sV7(s _ yfn 



1 - 



37r2 



12 Vmc^ 
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2„3 



27r^ 
45 



T 



(164) 



And in the case of weak interactions, when mc^ <S Tc, the fraction (163) has the form 

/ 2- X 3/2 

ris - 1 - ( — I - 



T J 



C(3/2) VT, 
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4i//3 



1 - 




^ 3/2" 



(165) 



The superfluid fraction Ug can be either larger or smaller than the condensate fraction no, 
depending on temperature, the strength of interactions, that is, on the gas parameter 7, 

and on the strength of disorder v. Increasing v leads to the simultaneous disappearance of 
the superfluid and condensate fractions at the same critical v^. through a flrst-order phase 
transition. This transition takes place between the superfluid phase, with no 7^ 0, 7^ 0, 
no 7^ 0, and un < 1, and the normal phase with no = 0, n^ = 0, uq — 0, and un — 1- 

15 Conclusion 

A self-consistent mean-flcld theory has been developed for Bose systems in random external 
potentials. The suggested approach makes it possible to consider arbitrarily strong inter- 
actions and an arbitrary strength of disorder. In general, the Bose system consists of the 
following components: the condensate fraction no, the normal fraction un-i the glassy frac- 
tion nc, and the superfluid fraction Ug- In the limit of asymptotically weak interactions 
and disorder, the known results are reproduced. When increasing the strength of disorder, 
a flrst-order phase transition occurs from the superfluid phase to the normal phase. For the 
class of models we considered we have found no pure Bose glass phase. The temperature 
for the occurrence of the first-order phase transition turns out to be lower than the critical 
temperature Tc of the second-order phase transition for a Bose system without disorder. The 
presence of disorder slightly lowers the sound velocity, but increases the density fiuctuations, 
the isothermal compressibility, and the structure factor. 

It is interesting that switching on disorder may lead to nonperturbative effects. For 
instance, the uniform ideal Bose gas is stochastically unstable with respect to infinitesimally 
small noise. Perturbation theory cannot be used to calculate the internal energy contributed 
by random fields. Nor is perturbation theory suflicicnt when analyzing the stability condition 
related to the minimization of the thermodynamic potential. 

The aim of present paper has been to develop an approach for considering Bose systems 
with any interaction strength and arbitrary strength of disorder and to describe the general 
properties of such systems. We have restricted ourselves to investigating those results that 
could be derived by analytic means. The overall quantitative study of the system properties 
requires to solve the intricate system of equations for functions of temperature T, density 
p, gas parameter 7, and noise parameter v. Such an investigation can be accomplished only 
numerically. In view of the length of the present paper, we prefer not to overload it further 
by these numerical calculations. They will be presented in separate publications. 
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Appendix A 

Throughout the paper, we have assumed that the external potential ^(r) possesses the 
Fourier transform for all k. The zero-momentum transform 



^0 = 1 C{r) dr 



exists in the sense that it is finite. This implies that the function ^(r) is integrable, such 
that 

j ^{r) dr < oo . 

For a function ^(r) integrable in that sense, the Riemann-Lebesgue theorem (for details see 
e.g. Refs. [57,58]) states that the related Fourier transforms is finite for all k and tends 
to zero as A; = |k| tends to infinity, < oo, 







{k oo) . 



How tends to zero depends on further properties of ^(r). For instance, if we can use that 
the function ,^(r) is finite at r = 0, that is. 



le(0)|<oo. 



then we have 



V 



< oo 



Replacing summation by integration in the standard manner, we get 



< oo . 



It follows from here that, when /c — > oo, then tends to zero at least as 

1 



(a > 3) . 



in three or more dimensions. The field is related to C,k by Eq. (75). Assuming that cpk 
has the same asymptotic behavior as ^fc, so that ipk tends to zero at large k —>■ oo not slower 
than 

1 

(a > 3) , 



k" 



we find solution (147), which confirms this assumption. With the field (fk possessing this 
asymptotic behavior and the form of ujk defined in Eq. (68), the ratio (pk/ {}^k + ^) tends to 
zero at large k not slower than 



cufc + A 



oc 



k' 



{n > 5) 



This justifies the transformation of Eq. (75) into Eq. (146). 
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Appendix B 



For the summation of scries (148), we used the method of self-similar factor approximants 
[54-56], which was shown to be more general and more accurate than the method of Fade 
approximants. In addition, the latter method is known to be not uniquely defined, in the 
sense that for each finite series of order k, there exists a table of Fade approximants P[m/n]: 
with M + N — k. Contrary to this, for each finite series of order k, there is just one 
factor approximant. This is why we prefer here to use the method of self-similar factor 
approximants. The construction of the factor approximants is done in the following way. 

Assume that for a function f{x), there is an asymptotic expansion at a; 0, so that 
f{x) can be represented by 

k 

fk{x) = fo{x) a„x" . 

n=0 

Since ao can always be included in fo{x), we may set, without the loss of generality, that 
Qo — 1. The self-similar factor approximant, extrapolating the finite series fk{x) to the 
region of arbitrary x, is defined as 

mx)^fo{x)i[{i+A.xr\ 

i=l 



where 

Nk 



k/2, k = 2,4,... 

{k + l)/2, ^ = 3,5,... 



The coefficients and exponents rii are given by the equations 

Y^mA^^Bn in^l,2,...,k) , 

i=l 

in which 

(n- 1)! ^-0 rfx" /o(a;) 

and Ai is set to one for odd k. As is evident, the parameters Ai, rii, and Bn, being calculated 
for each given k-th order approximation, depend on the approximation order k, so that 
Ai — Aik, rii — nik, and Bn — B^k- However, for the simplicity of notations, this order 
dependence is not shown explicitly. A detailed description of the method of self-similar 
factor approximants is given in Refs. [54-56]. 
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Appendix C 

Integral (160) can be regularized by means of the resummation regularization based 
on the self-similar approximation theory [51-53]. The procedure is as follows. One, first, 
introduces a cutoff L making the integral always converging, 

h{e)^j^ COth(£X Vl + x^) (^^^1^3/2 • 

Removing this cutoff would return us back to the integral 

I{e) = lim h{e) . 

i— »00 

Instead, we may represent the integral /l(£) as a series with the help of the expansion 



yielding the series 



with the coefficients 



coth^ = V -2'""' 
n=o (2^)! 



-1 oo 

^ n=0 



0'2n = 



•2n r ^ 
)! Jo (1 



(2n)! Jo (l + x2)2-n ' 
where are Bernoulli numbers. When L — > oo, then 

lim qq — — 

L— >oo 4 

and other coefficients behave as 

(2n)!(4n-l) ' 

with n = 1, 2, . . .. 

The series for can be represented with the help of the self-similar factor approxi- 

mants [54-56] as 

^ n=l 

where the parameters An and powers are uniquely defined from the re-expansion proce- 
dure, when 11(e) is expanded in powers of £^ and compared with the initial series for 
Then all An and are uniquely expressed through the coefficients a2n- For example, in 
lower orders, we have 

bl - 264 bl 
Ai = — ; , ai 



62 ' bl- 264 ' 

where 

^~ Oo ~ (2n)!(4n-l)7r 
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Then 

^ Qtt ' ^ 3157r 



etc. Because of this, 



1 4 280 1 



35 ' ^ ~ Qtt L ' 
and so on. 

Removing the cutoff in the self-similar approximant /£(£) by setting L — > oo, we use the 
limit 



lim (1 + 

>oo 



1/L" 



valid for any m > and n > 0. As a result, we have 

Substituting this into Eq. (159), we come to the same form of the dispersed heat (161) as 
obtained earlier. 

Another variant of the resummation regularization would be by summing the series for 
in the form of the self-similar exponential approximants [59,60]. This procedure gives 



Il{e) = ^ exp (he^ exp {h,e^ •••)), 



where again 6„ = On/ao- When setting L oo, we use the fact that the Bernoulli numbers 
are alternating in sign, so that 62, ^6; ^lO; and so on tend, polynomially in L, to plus infinity, 
while 64, 63, 612, and like that tend polynomially to minus infinity. Then we obtain 

' e ' 

which again leads to the same form (161). In this way, all considered variants of regularizing 
integral (160) give us the same expression (161), which confirms its general validity. 
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